Simulation of complete many-body quantum dynamics using controlled 

quantum— semiclassical hybrids 
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A controlled hybridization between full quantum dynamics and semiclassical approaches (mean- 
field and truncated Wigner) is implemented for interacting many-boson systems. It is then demon- 
strated how simulating the resulting hybrid evolution equations allows one to obtain the full quantum 
dynamics for much longer times than is possible using an exact treatment directly. A collision of 
sodium BECs with 1.5 x 10 J atoms is simulated, in a regime that is difficult to describe semi- 
classically. The uncertainty of physical quantities depends on the statistics of the full quantum 
prediction. Cutoffs are minimised to a discretization of the Hamiltonian. The technique presented 
is quite general and extension to other systems is considered. 



PACS numbers: 03.75.Kk, 05.30.-d, 05.10.Gg, 67.85.De 

The calculation of the full quantum dynamics of a 
many-body interacting system from the microscopic de- 
scription is a long-standing "difficult" problem with po- 
tential applications in many fields of physics — if only 
one could make it numerically tractable. The difficulty 
is that the size of the Hilbert space grows exponentially 
with the number of particles or orbitals, while path in- 
tegral Monte Carlo is foiled by the rapid appearance of 
random phases. How new headway against this problem 
can be made will be demonstrated below. 

Outside of fully integrable systems or ID, where 
MPS/DMRG-based methods are successful, simplified 
descriptions are used, e.g. mean-field theory, Bogoli- 
ubov diagonalization, long-wavelength or strong interac- 
tion expansions, and Wigner-distribution based "c-field" 
methods[l|, 0, 03]. However, some interesting problems 
fall outside the regimes of validity of these, typically 
where several competing effects are important or there 
is a transition between regimes that require different ap- 
proximations. In quantum gases this occurs with rising 
density when interactions between the coherent compo- 
nent and incoherent particles already become of essence 
during the evolution, but the gas is not yet dense enough 
for the c-field descriptions to describe it with only highly 
occupied modes. (See 0| for a comprehensive review of c- 
field methods and their validity). This may occur e.g. in 
quenches of the gas[3], colliding BECsjl, 0, 0], dynamics 
of the cooling and trapping, shock waves and the effects 
of obstacles[8| or disorder[9j. 

This kind of dynamics is often amenable to phase- 
space approaches that randomly sample the full 
quantum dynamics, such as positive-P[10|, stochastic 
wavefunctions 11| , and stochastic gauges [12j. They are 
successful when collective behaviour is important, but 
interactions between individual particles are not too 
strong. The density matrix p of the system is re-described 
in terms of a probability distribution p = J P(v)A(v)dv 
of basis operators A that is subsequently randomly sam- 



pled. These samples v are then evolved according to 
stochastic evolution equations that are chosen to keep 
the entire quantum dynamics of the microscopic descrip- 
tion. A serious limitation is the "noise catastrophe": Af- 
ter some finite time, an exponential (or faster) growth 
of the noise variance occurs, imposing a maximum fea- 
sible simulation time U \™ II 311 . While some phenomena 
can be simulated 1J, Uj£, UM , an extension of £ S im is much 



sought-after, and will be demonstrated here. 

The underlying reasons why phase-space methods can 
overcome the Hilbert space complexity, are that quanti- 
ties of physical interest usually involve contributions from 
many particles, and that limited precision is sufficient if 
it is well controlled. As in Monte-Carlo methods, there 
is no need to follow the amplitudes of all possible con- 
figurations as long as one can predict physical quantities 
with a well- controlled uncertainty. However — and now 
we come to the central idea to be demonstrated here — 
this can be taken further: There is also no true need to 
actually follow the troublesome exact quantum evolution 
equations provided that one can still predict what they 
would give with a well- controlled uncertainty. 

How can such a roundabout prediction be achieved? If 
one has at one's disposal two, or more, independent ap- 
proximate methods that produce evolution equations 11 A" 
and "23" without a noise catastrophe, but which bear suf- 
ficient resemblance to the full quantum dynamics equa- 
tions " Q", then hybrid equations can be constructed (pos- 
sibly ad-hoc) with a continuous blending parameter A in 
a scheme resembling 

H A = {1-\)A + \Q ; H B = (1 - A)23 + AQ. 

whose details will be non-universal. Here A = 1 gives full 
quantum dynamics, and A = the original approximate 
methods. The hybrids will still contain a noise catastro- 
phe, but at a later time than the full quantum treatment 
Q. Therefore, long times t > t^ m that are not accessible 
by Q will be accessible by some range of A € [0, A max (i) ]■ 
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If a physical quantity varies smoothly, preferably 
monotonically, as a function of A for hybrid 7iyi(A), then 
an extrapolation can be made to A = 1, based on sev- 
eral calculations in the accessible range [0, A max (i) < 1]. 
One extrapolation is not yet very convincing, however, 
it can be checked using the other independent hybrids 
7ie(A), .... When they all agree, one has an "interpo- 
lation between extrapolations" that is robust and much 
more reliable. Conceptually this step is similar to com- 
paring results obtained using different summation tech- 
niques in diagrammatic Monte-Carlo calculations [l7j. 

The remainder of this letter will demonstrate this pro- 
cedure on a system of colliding BECs (schematic shown 
in The parameters are chosen to be close to 

an early experiment at MIT0], but deliberately with 
fewer atoms, to put the system in the dilute yet Bose- 
stimulated regime where truncated Wigner and simple 
quasiparticle methods fail: An N = 1.5 x 10 5 atom BEC 
of 23 Na is prepared in an elongated magnetic trap with 
frequencies 20 x 80 x 80 Hz, at a temperature low enough 
to discount the thermal component (not unusual in ex- 
periments). A brief Bragg laser pulse coherently imparts 
a velocity kick of 2vq = 19.64mm/s to half the atoms 
along the long (x) condensate axis. The speed of the 
kicked atoms is supersonic (sound velocity in the cloud 
is < 3.1 mm/s). The trap is simultaneously turned off so 
that the wave-packets collide freely, producing a halo of 
scattered atom pairs moving at speeds ~ vq relative to 
the overall centre of mass. This scattered halo exhibits 
a rich behaviour, which has been the rep eated focus of 
experiments @, @, Q and theory [H, EE EE HE HJ H2, HI • 

The high-density regime of a similar system has been 
treated in detail with c-field methods in [23|]. Bogoliubov 
expansions and/or a pair-creation simplification treat the 
spontaneous regime, or special cases when BEC evolu 



tion is negligible or speed is highly supersonic|20l. l21 | 
(A stochastic Bogoliubov treatment gives promising re- 
sults in broader cases[24j]). However, major discrepan- 
cies between predictions for halo density and correla- 
tions arise when BEC evolution or Bose stimulation is 
appreciable. Correlations depend on the sizes of phase 
grains ji^l, which develop a complicated and poorly un- 
derstood shape HE HI and dynamics (13. EE H3] in this 
case. Parallels to unresolved questions in other fields of 
physics have been noted, such as the "HBT puzzle" in 
heavy ion collisions [25J. Trustworthy calculations that 
reach the end of the collision (observed in experiments^ 
but not reached by positive-P [IE E3|) could shed light 
on all these issues. 

Fig. [1] includes predictions from Gross-Pitaevskii (GP) 
mean field, truncated Wigner, and Positive-P calcula- 
tions. The time reachable by positive-P (t^ m ) is less than 
a half of the collision time i co i] ~ 1400/is, and both GP 
and Wigner give an error. The first does not treat scat- 
tering, while for a lattice fine enough to encompass all 
physics the second becomes valid only for N > 10 6 atoms 
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PIG. 1: Wigner (purple), positive-P (red), GP (dashed) and 
hybrid TLa calculations at various blending parameters A. (a): 
Total number of scattered atoms, from integration of k-space 
density (excluding the narrow BEC region), (b): Peak density 



of the halo (at v x 



0, v y = 9.37mm/s in velocity space) 



Triple lines show la uncertainty. 

(one needs > 0(1) atoms per lattice site[H). N.b. the 
/c-dependent difference between g and its effective lattice 
value[l| is < 3% here, so it has not been corrected for. 

Now let us turn to obtaining the full quantum dynam- 
ics for times longer than with the positive-P. The dynam- 
ics equations in the truncated Wigner, GP, and positive- 
P treatments share the GP kernel with certain additions, 
and turn out similar enough to play the role of the A, B, 
and Q. 

The dynamical GP equation for the complex field 
■0(x, t) corresponding to the cold atom Hamiltonian H = 
jd 3 x [$t( x )ff sp (x)$(x) + §§t( x )2$( x )2l ig ^( x ) = 

[H sp (x) + g\ip(x.)\ 2 ] An initial condensate wave- 

function 0gp(x) normalised to / d 3 x|^Gp(x)| 2 = N 
leads to initial conditions ^(x, 0) = </>gp(x). Expecta- 
tion values of observables (O) are calculated by making 
the replacements ^ — > ip an d ^ — ► ^* m O. For exam- 
ple, the density is n(x) = |-0(x)| 2 . 

In the truncated Wigner method, the dynamics is ob- 
tained by standard methods (e.g.[26|) based on the basis 
operator identities (x dependence implied) 
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whose importance for us will be seen below. The equa- 
tion of motion is as for GP but with the replacement 
\%jj\ 2 — > (|-0| 2 — 1) on the RHS. However, in the initial con- 
ditions the condensate field is admixed with half a virtual 
particle per mode as ^(x, 0) = ^gp(x) + ?7(x)/v / 2, where 
7?(x) is a local complex Gaussian noise with the ensemble 
averages (??(x)) = (r](x.)i](x.')) = and (r](x)7](x')*) = 
S 3 (x — x'). To calculate observables one ensemble av- 
erages a modified expression f[0] that is obtained via 



and subsequent re- 



(O) = Tr Op = j dvP(v) Tr OA 

placements ([1]), which give JdvP(v)f(v). E.g. n(x) = 
(l^(x)| 2 -i). 

The positive-P method uses two independent fields 
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t/>i(x, f) and ^(x, t) and the identities 



A$= 



^2 + 



V'i + off 



(2) 



The obey the Ito stochastic equations 



ifr^i(x) = [-Hsp(x) +gp(x) - 
ihip 2 (x) = [-ffsp(x) + gp(x)* 



V^Ci(x,i)] Vi(x) 
^^£2 (*,*)] ^2(x) 



(3) 



with "complex density" p(x) = -0i(x)^2(x)*. Here the 
£j are delta-correlated real Gaussian noise fields with the 
ensemble averages (£j(x, £)) = and (&(x, t)^- (x', t')) = 
5ijS(t — t')<5 3 (x — x'). Initial conditions are Vj( x , 0) = 
</>gp(x) and observables are obtained with the replace- 
ments \t — > ipi and "J^ — > ^2- 

The next step will be to hybridize the truncated 
Wigner with the positive-P into treatment Ha- It is 
most straightforward to proceed from hybrid operator 
identities for an off-diagonal expansion 
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(4) 



One obtains: n.(x) = (V'i(x) , 02(x)* — ^-o^) and initial 

ipj(x, 0) = 0gp(x) + T)(x)y~^^. The usual truncated- 
Wigner-like discarding of high-order derivatives in the 
relevant Fokker-Planck equations, gives dynamics 

ihipi(x) = [-ff S p(x) + c/p'(x) - v/^£i(x,*)] ^i(x) 
i^ 2 (x) = [-ffsp(x) +gp'(x)* - i^fig\&{-x,t)\ ^ 2 (x) 

with jo'(x) = p(x) + A— 1. As an aside, this corresponds to 
a representation based on an off-diagonal operator basis 
using s-ordered[i3| coherent-like states with s = A (See 
for details) . Fig. [TJ shows the performance of this 
hybrid for several values of A for two halo quantities of 
interest. As desired, A < 1 calculations last for longer 
than the full quantum dynamics. Here the simulation 
time scales as i S i m «oc 1/A, but this is not universal. 

Hybridization of the GP and positive-P methods into 
treatment Hb simply entails replacing ^/Tg by y/igX in 
the equations ([3]) and following the positive-P prescrip- 
tion from then on. Here t s - lm oc 1/A 2 . 

With hybrids in hand, extrapolations of the total num- 
ber of scattered atoms to the full QD limit A = 1 are 
shown in Fig. [2] for several times > Halo peak den- 
sity is in[l|3]- 

An issue here is deciding upon a fitting function - lin- 
ear, quadratic, otherwise? Firstly, an acceptable fit must 
not have any statistically significant mismatch with the 
data. Secondly, to exclude spurious ill-conditioned pa- 
rameters, one should choose a fit that minimises the un- 
certainty in the extrapolated value at A = 1 (see below). 
One must also beware of possible stiffness in the unseen 
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FIG. 2: A-dependent predictions for several times > t^ m 
(symbols) and corresponding quadratic fits (dashed line). Fit- 
ting is via minimisation of rms deviation in units of la data 
uncertainty. Data points use ~ 300 — 1000 trajectories. 

A, and sensitivity to this is the primary reason why sev- 
eral independent hybrids are needed. Details of Fig.[2]are 
consistent with a lack of stiffness in the unsimulated large 
A region: Firstly, for t at which the whole A sequence is 
seen, there are no inflections. Secondly, the two hybrids 
approach the A = 1 value from different sides but agree. 
Also, extrapolations from only a low-A portion of the 
available data should agree with ones that use the whole 
sequence. This is confirmed in (l8| . 

Agreement between the Ha and Hb extrapolations in 
Fig. [2] is rather good at long times, but it remains to 
provide a well-defined uncertainty for the final predic- 
tion. Methods to obtain the statistical uncertainty of 
the A = 1 extrapolation are known In this endeav- 
our it is very helpful to know the underlying distribution 
of the data points v(A), which are ensemble averaged 
observables. Conveniently, it is known to be Gaussian 
by the central limit theorem, and the shown la uncer- 
tainty Aw (A) is its standard deviation. One rather sim- 
ple way to proceed is to generate a number Ns ^S> 1 
of "synthetic" data sets, where in the jth set one gener- 
ates Vj(X) = v(X) + £j(A)Au(A), with £j being Gaussian 
random variables of variance 1, mean zero. The syn- 
thetic data Vj are distributed with the same mean as the 
original v but double the variance. Now one calculates 
an extrapolated QD prediction Vj(l) for A = 1 for each 
synthetic set j, and uses the distribution of these Vj(l) 
to obtain the final uncertainty Av(l). Predictions from 
Ha and Hb that match within statistical uncertainty are 
trustworthy to this accuracy. The final predictions from 
both hybrid methods for the number of scattered atoms 
are shown in Fig. 02, and for halo density in [lil ]. 

One sees that the useful simulation time has been ex- 
tended several-fold, allows one to reach the end of the 
collision here, and determine the total scattered atoms 
to be 8800 ± 400 (at t= 1.7ms). The much worse preci- 
sion of the Ha result stems from the inherent vacuum 
noise in Wigner calculations and shorter segment of A 
values. However, for halo density, it is Hb that is more 
noisy. 

Regarding limits of applicability, at very long times the 
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FIG. 3: Predictions of from hybrids Ha and Hb compared 
with short-time full quantum dynamics and approximate 
methods. Triple lines, where visible, are lcr uncertainty. Uses 
« 10 — 20 values of A, as per Fig. [2] 

uncertainty becomes excessive for all hybrids since the 
short A intervals give badly conditioned extrapolations. 
Hence the bare simulation time in the Q treatment must 
not be too small to ensure a sufficiently long A interval. 
It is also crucial that the blending A enter the dynam- 
ics in a global way: Artificial boundaries^, could 
make observables depend stiffly on the boundary posi- 
tion. For cold gases low densities can be treated pertur- 
batively, while at high enough densities c-field treatments 
are valid, so that one expects that the blending method 
will be most useful at intermediate densities that "fall 
through the cracks" between these two methods. The 
relative simplicity of not requiring a projection onto low- 
energy modes may also make blending appealing in other 
regimes. 

Finally, while the emphasis has been on cold bo- 
son dynamics, the general equation-blending approach 
should be broadly applicable. For hard-core boson or 
fermion systems other approximations would have to be 
hybridised with a different complete phase-space descrip- 
tion Q. One can also hybridise "imaginary-time" evolu- 
tion for thermal equilibrium states, or Monte-Carlo path- 
integrals with the aim of predicting the ab-initio result for 
longer j3 = l/T than is normally allowed by the fermion 
sign problem. 

Concluding, it has been demonstrated how the full 
quantum dynamics of a macroscopic interacting 3D sys- 
tem can be calculated for much longer times than was 
possible with the previously most effective method, the 
positive-P representation. Quantitative predictions for 
BEC collisions in the dilute stimulated regime were ob- 
tained. The hybrid dynamical equations used, while not 
actually simulating complete quantum dynamics per se, 
can be used to confidently predict the full quantum dy- 
namics (within a given accuracy) when several families 
of hybrids are available. 
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The BEC collision 
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The system simulated, (a): Schematic of the BEC collision 
in real space in the lab frame, (b): Slice of the velocity 
distribution p in the center-of-mass frame at v z = and 
t = 670/is calculated using the positive-P method. This is 
about a third of the collision time, and the maximum time 
achievable with that method. The condensates are located 
around v x = ±vq = ±9.82mm/s. The halo of scattered atoms 
is clearly seen, as are the coherent frequency doubling peaks 
at ±3uq ~ 30mm/s. The collision is along the x axis. 



The relationship of the hybrid Ha to s-ordered 
operators 

First, a brief exposition of the standard formalism used 
in deriving phase-space quantum dynamics will be neces- 
sary. Writing the state of the system as a density matrix 
p, it can also be expressed as a distribution 



p= / dvP{v)A(v). 



(5) 



over a family of basis operators A(v) parameterised by 
variables in the set v. If the distribution P(v) is real and 
non-negative, this corresponds, in turn, to an ensemble 
of iS sets of random variables v ("configurations") chosen 
according to the distribution P, in the limit when S — > 
oo. In practice one computes a finite but large ensemble 
(S ^> 1) and knows properties of p to within a statistical 
uncertainty that can be confidently estimated from the 
properties of the finite ensemble. 

The dynamics of the system is described by the master 
equation 



lTi dt 



H,p 



while expectation values of observables are 



(O) = Tr Op 



(6) 



(7) 



These are most readily related to the computational en- 
semble of random variables through the use of the "oper- 
ator identities", that are specific to each formulation. 

For example, in the positive-P method one chooses A 
to be an off-diagonal coherent-state operator. Letting 
x label discrete points in the computational lattice with 
AV volume per point, defining 



one has 



App(<7) = n 

X 

where v — {ai, 0.2}^ 



a j (y)=il> j {x)/VAV, 

ai(x)) x (a 2 (x)|, 



(a 2 (x)| x |ai(x))> 



(8) 



M 2 /2 «at 



10), 



is a coherent state on the x lattice point with the complex 
amplitude a and anihilation operator a x = 4'(x)\/ AV. 
Then, one finds (omitting ubiquitous local x dependence) 
the operator identities: 



i&App = i'lApp 



pp 



A PP ^ = 



d_ 

dtpi 
d 
W2 



V>2 + — 

i'l- 



App 

App, 



which are the source of the positive-P identities in the 
main text. Combined with §5§ and (J6J these allow one 
to obtain a partial differential equation for P(v, t) that is 
equivalent to the full quantum evolution of p{i). For the 
positive-P representation, this is a Fokker-Planck equa- 
tion, and it corresponds exactly to the Langevin equa- 
tions given in (5) of the main text Combining the iden- 



tities with (7]) and Tr 



(O) 



A 



pp 



= 1 one finds 



P(v)fo(v)dv 



with a function fo that is obtained from O via the oper- 
ator identities, so that in the calculation it corresponds 
to an ensemble average of fo- For example, for O = 
$t( x )$(x), the function is 1 f = '0|(x)Vi(x). The ini- 
tial coherent state corresponds to P = Y[ x j S 1 - 3 ^ (tpj(x) — 
0Gp(x)). 

It has been shown 2 that the Glauber-Sudarshan P dis- 
tribution described by a coherent state operator basis 



A G spw=nH x )>xMx) 



(similar to the positive-P but diagonal) can be described 
as the limit of a representation over s-ordered basis states 

Agsp = lim A s 

s^l- 

where s can take on continuous values from -1 to 1, and 



Tr[5(a) x f(0,- S ) x 5-i(«)> 



Here 



v ' ' 1 + s V 1 + s 



is a kernel operator that becomes the vacuum |0)(0| in 
the limit of s — > 1~ and the local displacement operator 
is 



D(a) 



so that coherent states are |a) = D(a)|0). It was also 
shown there that the Wigner distribution corresponds to 
s = 0, hence a variation of s from to 1 looks like a good 
candidate to create the Jij^, hybrid formulation between 
truncated Wigner and positive-P. The "truncation" refers 



1 So = ipi {~ x )4>2 (x) can also be obtained, but gives the same value 
of (O) in the S — > oo limit. 

2 K. E. Cahill and R. J. Glauber, Phys. Rev. 177, 1857 (1969); 
ibid. 177, 1882 (1969) 
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to ad-hoc removal of third order 3 partial derivatives of 
the Wigner distribution P in its evolution equation to 
make it interpretable as Langevin stochastic equations of 
the samples. This removal is the reason why truncated 
Wigner treatments do not reproduce the full quantum 
dynamics. 

First, though, one must take into account the off- 
diagonality that is responsible for the difference between 
the Glauber-Sudarshan P and positive-P: A PP ^ Aqs P . 
Notably one of the bases 4 that reproduces the positive-P 
is 



One obtains the identities 5 



A 



pp 



>i = n 



d(v)xT(0,-l) x d-V) 



(10) 



Tr [d(fl)xT(0,-l)xd- 1 (d)> 
= \{d{v) x f ({),-!) J- 1 {v)^ 

X 

where the "displacement-like" operator 



is obtained by the replacement a — ► a% , a* — ► a| in D(a), 
and the second line follows because the trace in the de- 
nominator evaluates to one. The reason for this partic- 
ular replacement is that for the positive-P distribution 
one requires A to depend analytically on two separate 
complex variables, hence their complex conjugates must 
be removed. Here these analytic variables are <x\ and ot\. 

The extension of this A onto a family of s-ordered bases 

is 



&A A 
A?& 
A A ty 



1>i 

02 
01 



1 - 


s d 


2 


dtp* 


1 + 


s d 


2 




1 - 


s d 


2 




1 + 


s d 


2 


or 2 



Af 

A A 
A A 



which are exactly the same as was obtained by a naive 
blending of the operator identities in the main text pro- 
vided we identify A = s. 

Regarding initial conditions, the diagonal s-ordered 
representation © for a coherent state \4>gp) was found 
by Cahill and Glauber to be Gaussian 



p w = IL — cx p 

XJ - 1 — s 



2|^(x)-cfcp(x)p 
AV(1 - s) 



(12) 



When one additionally imposes ipi = -02 = "0 as is done 
in the main text, this is equivalent to ifTTj) . justifying 
the initial conditions given in the main text that contain 
complex Gaussian noise of variance (1 — s)/2. 



a>i = n 



Tr [2(^(0,-8) Jh^if) 



(11) 



This then interpolates towards the Wigner representa- 
tion. Note that since the truncated Wigner evolution is 
deterministic, then if one takes the formally off-diagonal 
basis set with s = but imposes 5(ipi — 1P2) in the ini- 
tial conditions, it will remain exactly equivalent to the 
normal truncated Wigner formulation of (JH) with s = 0. 



3 And higher order terms if necessary, although for the cold atom 
Hamiltonian considered in this letter, only partial derivatives up 
to third order are present in the Wigner representation. 

4 Though not the only one. jDther^ ways of writing A such as e.g. 
B(ai)f(0, -l)D(a*)/Tt[D( ai )f{0,-l)D(a*)] can also repro- 
duce the positive-P formulation but are not useful for general- 
isation to s < 1, and do not reproduce the same itermediate 
operator identities. 



5 For example, by comparison of expressions for LHS and RHS 
when T(0, — s) is expanded in number states. 
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Halo density calculations 



Extrapolation from partial A segment 



t=1010ns 

t=1895|is 
. 1 -5 h circles: H :tWigner/+P 
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JU- ' 
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A-dependent predictions of halo density (at v x = v z = 0, 
u H = 9.37mm/s in velocity space) for several times (circles) 
with uncertainty shown as vertical bars at the same location. 
The corresponding fits (dashed) are quadratic for the TLs hy- 
brid, and constant- value for TLa- Fitting is via minimisation 
of rms deviation in units of la data uncertainty. Linear or 
quadratic fits to the TLa hybrid data are not more statisti- 
cally significant than the constant- value fit, and hence would 
be poorly conditioned. 
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Predictions of the number of scattered atoms at several times, 
as a function of the A segment A 6 [0, A max ] used for extrapo- 
lation from a quadratic fit to Ha results. Triple lines, where 
visible, are la uncertainty. Dashed lines indicate the final pre- 
dictions using all the available A values. Data used was from 
the same simulations as in Fig. 2 of the main text. There is 
no statistically significant trend with A max visible, suggesting 
that the fitting function that is a quadratic polynomial in A 
is appropriate within statistical precision. 





t [units of ms] 
Predictions of halo density (at v x 



0, 



9.37mm/s 



in velocity space) from hybrids Ha and Tie compared with 
short-time full quantum dynamics and approximate methods. 
Triple lines, where visible, are la uncertainty. Prediction data 
based on w 10 - 20 values of A, each with « 300 - 1000 
trajectories, and quadratic / constant-value fitting for TLa / 
TLb hybrids, respectively. Note the agreement with truncated 
Wigner to within statistical uncertainty. Times detailed in 
the previous figure (above) are highlighted. 



